Water Rights Restored to the Gila River¶

The impacts of irrigation on vegetation health in the Gila River Valley¶

Background¶

Observing vegetation health from space¶

We will look at vegetation health using NDVI (Normalized Difference Vegetation Index). How does it work? First, we need to learn about spectral reflectance signatures.

Every object reflects some wavelengths of light more or less than others. We can see this with our eyes, since, for example, plants reflect a lot of green in the summer, and then as that green diminishes in the fall they look more yellow or orange. The image below shows spectral signatures for water, soil, and vegetation:

> Image source: SEOS Project

Healthy vegetation reflects a lot of Near-InfraRed (NIR) radiation. Less healthy vegetation reflects a similar amounts of the visible light spectra, but less NIR radiation. We don’t see a huge drop in Green radiation until the plant is very stressed or dead. That means that NIR allows us to get ahead of what we can see with our eyes.

Different species of plants reflect different spectral signatures, but the pattern of the signatures across species and sitations is similar. NDVI compares the amount of NIR reflectance to the amount of Red reflectance, thus accounting for many of the species differences and isolating the health of the plant. The formula for calculating NDVI is:

$$NDVI = \frac{(NIR - Red)}{(NIR + Red)}$$

Read More

Read more about NDVI and other vegetation indices:

  • earthdatascience.org
  • USGS
In [1]:
id = 'stars'
site_name = 'Gila River Indian Community'
project_name = 'Gila River Vegetation'
boundary_dir = 'tl_2020_us_aitsn'
event = 'water rights case'
start_year = '2001'
end_year = '2022'
event_year = '2012'

00: Setup¶

Import Libraries and Download Sample Data¶

In [2]:
import json
import os
import pathlib

from glob import glob #C Combine data arrays

import earthpy # Manage local data
import geopandas as gpd # Work with geospatial vector data
import pandas as pd # Work with tabular data
import hvplot.pandas


import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point

import hvplot.xarray
import rioxarray as rxr # To work with raster data
import xarray as xr # To work with data arrays and raster data
In [3]:
project = earthpy.Project(
    'Gila River Vegetation', dirname='Vegetation')
project.get_data()

01: Create Site Map¶

Study Area: Gila River Indian Community¶

Earth Data Science data formats¶

In Earth Data Science, we get data in three main formats:

Data type Descriptions Common file formats Python type
Time Series The same data points (e.g. streamflow) collected multiple times over time Tabular formats (e.g. .csv, or .xlsx) pandas DataFrame
Vector Points, lines, and areas (with coordinates) Shapefile (often an archive like a .zip file because a Shapefile is actually a collection of at least 3 files) geopandas GeoDataFrame
Raster Evenly spaced spatial grid (with coordinates) GeoTIFF (.tif), NetCDF (.nc), HDF (.hdf) rioxarray DataArray
Read More

about vector data and raster data in the textbook.

Reflect and Respond¶

  • For this coding challenge, we are interested in the Gila River Indian Community (GRIC) boundary and the health of the vegetation in the area measure on a scale from -1 to 1.

What type of data do you think the boundary will be?

  • The data for the boundary will probably be a polygon data type.

What type of data do you think the vegetation health will be?

  • The vegetation health will probably be an object (variable) data type.

Load the Gila River Indian Community boundary¶

In [4]:
# Load in the boundary data
aitsn_gdf = gpd.read_file(project.project_dir / 'tl_2020_us_aitsn')

# Check that it worked
aitsn_gdf.head()
Out[4]:
AIANNHCE TRSUBCE TRSUBNS GEOID NAME NAMELSAD LSAD CLASSFP MTFCC FUNCSTAT ALAND AWATER INTPTLAT INTPTLON geometry
0 2430 653 02419073 2430653 Red Valley Red Valley Chapter T2 D7 G2300 A 922036695 195247 +36.6294607 -109.0550394 POLYGON ((-109.2827 36.64644, -109.28181 36.65...
1 2430 665 02419077 2430665 Rock Point Rock Point Chapter T2 D7 G2300 A 720360268 88806 +36.6598701 -109.6166836 POLYGON ((-109.85922 36.49859, -109.85521 36.5...
2 2430 675 02419081 2430675 Rough Rock Rough Rock Chapter T2 D7 G2300 A 364475668 216144 +36.3976971 -109.7695183 POLYGON ((-109.93053 36.40672, -109.92923 36.4...
3 2430 325 02418975 2430325 Indian Wells Indian Wells Chapter T2 D7 G2300 A 717835323 133795 +35.3248534 -110.0855000 POLYGON ((-110.24222 35.36327, -110.24215 35.3...
4 2430 355 02418983 2430355 Kayenta Kayenta Chapter T2 D7 G2300 A 1419241065 1982848 +36.6884391 -110.3045616 POLYGON ((-110.56817 36.73489, -110.56603 36.7...
In [5]:
# Make interactive plot to view data and find site location

aitsn_gdf.hvplot(
    geo=True, tiles='EsriImagery', 
    frame_width=500,
    legend=False, fill_color=None, edge_color='white',
    # This parameter makes all the column values in the dataset visible.
    hover_cols='all')
WARNING:param.main: edge_color option not found for polygons plot with bokeh; similar options include: ['muted_color', 'bgcolor', 'line_color']
Out[5]:

Reflect and Respond¶

What column could you use to uniquely identify the subdivisions of the reservation you want to study using this interactive map? What value do you need to use to filter the GeoDataFrame?¶

  • The column AIANNHCE is the identifier for the tribal subdivision
  • Looking for Phoenix, Arizona at lat,lon 33 N, 112 W in the interactive map, we find AIANNHCE is number 1310
In [6]:
# Check what type of data the AIANNHCE column is
aitsn_gdf.info()
<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 484 entries, 0 to 483
Data columns (total 15 columns):
 #   Column    Non-Null Count  Dtype   
---  ------    --------------  -----   
 0   AIANNHCE  484 non-null    object  
 1   TRSUBCE   484 non-null    object  
 2   TRSUBNS   484 non-null    object  
 3   GEOID     484 non-null    object  
 4   NAME      484 non-null    object  
 5   NAMELSAD  484 non-null    object  
 6   LSAD      484 non-null    object  
 7   CLASSFP   484 non-null    object  
 8   MTFCC     484 non-null    object  
 9   FUNCSTAT  484 non-null    object  
 10  ALAND     484 non-null    int64   
 11  AWATER    484 non-null    int64   
 12  INTPTLAT  484 non-null    object  
 13  INTPTLON  484 non-null    object  
 14  geometry  484 non-null    geometry
dtypes: geometry(1), int64(2), object(12)
memory usage: 56.8+ KB

Reflect and Respond¶

What is the data type of the AIANNHCE column?¶

How will that affect your code?

  • The AIANNHCE column is a geodata frame of objects. This means that we will have to use geopandas to work with a dataframe array.
In [7]:
# Select and merge the Gila River subdivisions

gila_gdf = aitsn_gdf.loc[aitsn_gdf.AIANNHCE=='1310'].dissolve()

# Plot the results with web tile images

gila_gdf.hvplot(
    geo=True, tiles='EsriImagery',
    fill_color=None, line_color='black',
    title='Gila River Indian Community',
    frame_width=500)
Out[7]:
In [8]:
%store project gila_gdf aitsn_gdf
Stored 'project' (Project)
Stored 'gila_gdf' (GeoDataFrame)
Stored 'aitsn_gdf' (GeoDataFrame)

02: Wrangle Data¶

Load in NDVI data¶

  • In this section, we need to extract the date from the filename. We will use the glob.rglob() function in Python to find all files that match a pattern.

Reflect and Respond¶

Take a look at the file names for the NDVI files. What do you notice is the same for all the files?

Keep in mind that you only want the NDVI files, not the Quality files (which will flag potential incorrect measurements)

  • Looking in the gila-river-ndvi folder, we can see that all the file names all include NDVI and are of type .tif, so we can use the wildcard character (), for files 'NDVI*.tif'
In [9]:
# Get a sorted list of NDVI tif file paths
ndvi_paths = sorted(list(project.project_dir.rglob('*NDVI*.tif')))

# See how many files are in this list
print (len(ndvi_paths)) 

# Display the first and last three files paths to check the pattern
(ndvi_paths[:3], ndvi_paths[-3:])
154
Out[9]:
([PosixPath('/workspaces/data/vegetation/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001145000000_aid0001.tif'),
  PosixPath('/workspaces/data/vegetation/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001161000000_aid0001.tif'),
  PosixPath('/workspaces/data/vegetation/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001177000000_aid0001.tif')],
 [PosixPath('/workspaces/data/vegetation/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022209000000_aid0001.tif'),
  PosixPath('/workspaces/data/vegetation/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022225000000_aid0001.tif'),
  PosixPath('/workspaces/data/vegetation/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022241000000_aid0001.tif')])
In [10]:
# Extract information from the data sets
doy_start = -25
doy_end = -19

# Loop through each NDVI image

ndvi_das = [] # Create an empty list of data arrays

for ndvi_path in ndvi_paths: # ndvi_path equivalent to i
    
    # Get date from file name
    doy = ndvi_path.name[doy_start:doy_end]
    date = pd.to_datetime(doy, format='%Y%j') # Use pandas to change doy to date
    
    # Open dataset
    da = rxr.open_rasterio(ndvi_path, mask_and_scale=True).squeeze() # Create data array dataset

    # Add date dimension and clean up metadata
    da = da.assign_coords({'date': date})
    da = da.expand_dims({'date': 1})
    da.name = 'NDVI'

    # Prepare for concatenation
    ndvi_das.append(da) # Add the data array in each loop to the list of NDVI arrays

# Return length of list
len(ndvi_das)
Out[10]:
154
In [11]:
# Combine NDVI images (Raster data) from all dates
ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
ndvi_da
/tmp/ipykernel_2734/570708293.py:2: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly.
  ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
/tmp/ipykernel_2734/570708293.py:2: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly.
  ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
Out[11]:
<xarray.Dataset> Size: 48MB
Dimensions:      (date: 154, y: 203, x: 382)
Coordinates:
    band         int64 8B 1
  * x            (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
  * y            (y) float64 2kB 33.39 33.39 33.38 33.38 ... 32.97 32.97 32.97
    spatial_ref  int64 8B 0
  * date         (date) datetime64[ns] 1kB 2001-01-14 2001-01-16 ... 2022-01-24
Data variables:
    NDVI         (date, y, x) float32 48MB 0.8282 0.6146 ... 0.2146 0.2085
xarray.Dataset
    • date: 154
    • y: 203
    • x: 382
    • band
      ()
      int64
      1
      array(1)
    • x
      (x)
      float64
      -112.3 -112.3 ... -111.5 -111.5
      array([-112.309375, -112.307292, -112.305208, ..., -111.519792, -111.517708,
             -111.515625], shape=(382,))
    • y
      (y)
      float64
      33.39 33.39 33.38 ... 32.97 32.97
      array([33.388542, 33.386458, 33.384375, ..., 32.971875, 32.969792, 32.967708],
            shape=(203,))
    • spatial_ref
      ()
      int64
      0
      crs_wkt :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      inverse_flattening :
      298.257223563
      reference_ellipsoid_name :
      WGS 84
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      geographic_crs_name :
      WGS 84
      horizontal_datum_name :
      World Geodetic System 1984
      grid_mapping_name :
      latitude_longitude
      spatial_ref :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      GeoTransform :
      -112.31041665660531 0.0020833333331466974 0.0 33.38958333034212 0.0 -0.0020833333331466974
      array(0)
    • date
      (date)
      datetime64[ns]
      2001-01-14 ... 2022-01-24
      array(['2001-01-14T00:00:00.000000000', '2001-01-16T00:00:00.000000000',
             '2001-01-17T00:00:00.000000000', '2001-01-19T00:00:00.000000000',
             '2001-01-20T00:00:00.000000000', '2001-01-22T00:00:00.000000000',
             '2001-01-24T00:00:00.000000000', '2002-01-14T00:00:00.000000000',
             '2002-01-16T00:00:00.000000000', '2002-01-17T00:00:00.000000000',
             '2002-01-19T00:00:00.000000000', '2002-01-20T00:00:00.000000000',
             '2002-01-22T00:00:00.000000000', '2002-01-24T00:00:00.000000000',
             '2003-01-14T00:00:00.000000000', '2003-01-16T00:00:00.000000000',
             '2003-01-17T00:00:00.000000000', '2003-01-19T00:00:00.000000000',
             '2003-01-20T00:00:00.000000000', '2003-01-22T00:00:00.000000000',
             '2003-01-24T00:00:00.000000000', '2004-01-14T00:00:00.000000000',
             '2004-01-16T00:00:00.000000000', '2004-01-17T00:00:00.000000000',
             '2004-01-19T00:00:00.000000000', '2004-01-20T00:00:00.000000000',
             '2004-01-22T00:00:00.000000000', '2004-01-24T00:00:00.000000000',
             '2005-01-14T00:00:00.000000000', '2005-01-16T00:00:00.000000000',
             '2005-01-17T00:00:00.000000000', '2005-01-19T00:00:00.000000000',
             '2005-01-20T00:00:00.000000000', '2005-01-22T00:00:00.000000000',
             '2005-01-24T00:00:00.000000000', '2006-01-14T00:00:00.000000000',
             '2006-01-16T00:00:00.000000000', '2006-01-17T00:00:00.000000000',
             '2006-01-19T00:00:00.000000000', '2006-01-20T00:00:00.000000000',
             '2006-01-22T00:00:00.000000000', '2006-01-24T00:00:00.000000000',
             '2007-01-14T00:00:00.000000000', '2007-01-16T00:00:00.000000000',
             '2007-01-17T00:00:00.000000000', '2007-01-19T00:00:00.000000000',
             '2007-01-20T00:00:00.000000000', '2007-01-22T00:00:00.000000000',
             '2007-01-24T00:00:00.000000000', '2008-01-14T00:00:00.000000000',
             '2008-01-16T00:00:00.000000000', '2008-01-17T00:00:00.000000000',
             '2008-01-19T00:00:00.000000000', '2008-01-20T00:00:00.000000000',
             '2008-01-22T00:00:00.000000000', '2008-01-24T00:00:00.000000000',
             '2009-01-14T00:00:00.000000000', '2009-01-16T00:00:00.000000000',
             '2009-01-17T00:00:00.000000000', '2009-01-19T00:00:00.000000000',
             '2009-01-20T00:00:00.000000000', '2009-01-22T00:00:00.000000000',
             '2009-01-24T00:00:00.000000000', '2010-01-14T00:00:00.000000000',
             '2010-01-16T00:00:00.000000000', '2010-01-17T00:00:00.000000000',
             '2010-01-19T00:00:00.000000000', '2010-01-20T00:00:00.000000000',
             '2010-01-22T00:00:00.000000000', '2010-01-24T00:00:00.000000000',
             '2011-01-14T00:00:00.000000000', '2011-01-16T00:00:00.000000000',
             '2011-01-17T00:00:00.000000000', '2011-01-19T00:00:00.000000000',
             '2011-01-20T00:00:00.000000000', '2011-01-22T00:00:00.000000000',
             '2011-01-24T00:00:00.000000000', '2012-01-14T00:00:00.000000000',
             '2012-01-16T00:00:00.000000000', '2012-01-17T00:00:00.000000000',
             '2012-01-19T00:00:00.000000000', '2012-01-20T00:00:00.000000000',
             '2012-01-22T00:00:00.000000000', '2012-01-24T00:00:00.000000000',
             '2013-01-14T00:00:00.000000000', '2013-01-16T00:00:00.000000000',
             '2013-01-17T00:00:00.000000000', '2013-01-19T00:00:00.000000000',
             '2013-01-20T00:00:00.000000000', '2013-01-22T00:00:00.000000000',
             '2013-01-24T00:00:00.000000000', '2014-01-14T00:00:00.000000000',
             '2014-01-16T00:00:00.000000000', '2014-01-17T00:00:00.000000000',
             '2014-01-19T00:00:00.000000000', '2014-01-20T00:00:00.000000000',
             '2014-01-22T00:00:00.000000000', '2014-01-24T00:00:00.000000000',
             '2015-01-14T00:00:00.000000000', '2015-01-16T00:00:00.000000000',
             '2015-01-17T00:00:00.000000000', '2015-01-19T00:00:00.000000000',
             '2015-01-20T00:00:00.000000000', '2015-01-22T00:00:00.000000000',
             '2015-01-24T00:00:00.000000000', '2016-01-14T00:00:00.000000000',
             '2016-01-16T00:00:00.000000000', '2016-01-17T00:00:00.000000000',
             '2016-01-19T00:00:00.000000000', '2016-01-20T00:00:00.000000000',
             '2016-01-22T00:00:00.000000000', '2016-01-24T00:00:00.000000000',
             '2017-01-14T00:00:00.000000000', '2017-01-16T00:00:00.000000000',
             '2017-01-17T00:00:00.000000000', '2017-01-19T00:00:00.000000000',
             '2017-01-20T00:00:00.000000000', '2017-01-22T00:00:00.000000000',
             '2017-01-24T00:00:00.000000000', '2018-01-14T00:00:00.000000000',
             '2018-01-16T00:00:00.000000000', '2018-01-17T00:00:00.000000000',
             '2018-01-19T00:00:00.000000000', '2018-01-20T00:00:00.000000000',
             '2018-01-22T00:00:00.000000000', '2018-01-24T00:00:00.000000000',
             '2019-01-14T00:00:00.000000000', '2019-01-16T00:00:00.000000000',
             '2019-01-17T00:00:00.000000000', '2019-01-19T00:00:00.000000000',
             '2019-01-20T00:00:00.000000000', '2019-01-22T00:00:00.000000000',
             '2019-01-24T00:00:00.000000000', '2020-01-14T00:00:00.000000000',
             '2020-01-16T00:00:00.000000000', '2020-01-17T00:00:00.000000000',
             '2020-01-19T00:00:00.000000000', '2020-01-20T00:00:00.000000000',
             '2020-01-22T00:00:00.000000000', '2020-01-24T00:00:00.000000000',
             '2021-01-14T00:00:00.000000000', '2021-01-16T00:00:00.000000000',
             '2021-01-17T00:00:00.000000000', '2021-01-19T00:00:00.000000000',
             '2021-01-20T00:00:00.000000000', '2021-01-22T00:00:00.000000000',
             '2021-01-24T00:00:00.000000000', '2022-01-14T00:00:00.000000000',
             '2022-01-16T00:00:00.000000000', '2022-01-17T00:00:00.000000000',
             '2022-01-19T00:00:00.000000000', '2022-01-20T00:00:00.000000000',
             '2022-01-22T00:00:00.000000000', '2022-01-24T00:00:00.000000000'],
            dtype='datetime64[ns]')
    • NDVI
      (date, y, x)
      float32
      0.8282 0.6146 ... 0.2146 0.2085
      units :
      NDVI
      AREA_OR_POINT :
      Area
      array([[[0.8282    , 0.6146    , 0.3796    , ..., 0.1542    ,
               0.1542    , 0.1774    ],
              [0.52349997, 0.32029998, 0.46359998, ..., 0.145     ,
               0.145     , 0.145     ],
              [0.52349997, 0.38509998, 0.4567    , ..., 0.0962    ,
               0.058     , 0.0641    ],
              ...,
              [0.13759999, 0.13759999, 0.1283    , ..., 0.1489    ,
               0.1489    , 0.1568    ],
              [0.1299    , 0.1299    , 0.12379999, ..., 0.1831    ,
               0.1831    , 0.23169999],
              [0.1347    , 0.1309    , 0.1329    , ..., 0.15709999,
               0.2384    , 0.2452    ]],
      
             [[0.5146    , 0.42249998, 0.3706    , ..., 0.1656    ,
               0.1656    , 0.1656    ],
              [0.3543    , 0.4138    , 0.3706    , ..., 0.1275    ,
               0.1275    , 0.0968    ],
              [0.48029998, 0.3543    , 0.3687    , ..., 0.075     ,
               0.075     , 0.0968    ],
      ...
              [0.1303    , 0.1303    , 0.1303    , ..., 0.1739    ,
               0.1739    , 0.1859    ],
              [0.13319999, 0.13319999, 0.1303    , ..., 0.19399999,
               0.19399999, 0.2383    ],
              [0.1662    , 0.13319999, 0.1348    , ..., 0.17289999,
               0.2034    , 0.23249999]],
      
             [[0.5902    , 0.3793    , 0.3145    , ..., 0.15359999,
               0.15359999, 0.1564    ],
              [0.6694    , 0.4878    , 0.32909998, ..., 0.159     ,
               0.159     , 0.1331    ],
              [0.3103    , 0.39409998, 0.4527    , ..., 0.154     ,
               0.1026    , 0.0448    ],
              ...,
              [0.1338    , 0.1338    , 0.1192    , ..., 0.174     ,
               0.174     , 0.163     ],
              [0.12889999, 0.12889999, 0.1205    , ..., 0.2146    ,
               0.2146    , 0.2025    ],
              [0.13069999, 0.1157    , 0.12629999, ..., 0.17359999,
               0.2146    , 0.2085    ]]], shape=(154, 203, 382), dtype=float32)
    • x
      PandasIndex
      PandasIndex(Index([-112.30937498993873, -112.30729165660559, -112.30520832327244,
             -112.30312498993929, -112.30104165660615,   -112.298958323273,
             -112.29687498993985,  -112.2947916566067, -112.29270832327356,
             -112.29062498994041,
             ...
             -111.53437499000816, -111.53229165667501, -111.53020832334187,
             -111.52812499000872, -111.52604165667557, -111.52395832334243,
             -111.52187499000928, -111.51979165667613, -111.51770832334299,
             -111.51562499000984],
            dtype='float64', name='x', length=382))
    • y
      PandasIndex
      PandasIndex(Index([ 33.38854166367555,   33.3864583303424,  33.38437499700925,
              33.38229166367611,  33.38020833034296,  33.37812499700981,
             33.376041663676666,  33.37395833034352,  33.37187499701037,
             33.369791663677226,
             ...
             32.986458330378234,  32.98437499704509,  32.98229166371194,
             32.980208330378794,  32.97812499704565,   32.9760416637125,
             32.973958330379354,  32.97187499704621,  32.96979166371306,
             32.967708330379914],
            dtype='float64', name='y', length=203))
    • date
      PandasIndex
      PandasIndex(DatetimeIndex(['2001-01-14', '2001-01-16', '2001-01-17', '2001-01-19',
                     '2001-01-20', '2001-01-22', '2001-01-24', '2002-01-14',
                     '2002-01-16', '2002-01-17',
                     ...
                     '2021-01-20', '2021-01-22', '2021-01-24', '2022-01-14',
                     '2022-01-16', '2022-01-17', '2022-01-19', '2022-01-20',
                     '2022-01-22', '2022-01-24'],
                    dtype='datetime64[ns]', name='date', length=154, freq=None))
In [12]:
%store ndvi_da ndvi_das ndvi_paths
Stored 'ndvi_da' (Dataset)
Stored 'ndvi_das' (list)
Stored 'ndvi_paths' (list)

03: Plot Data¶

In [13]:
# Plot first .tif image
rxr.open_rasterio(ndvi_paths[0], mask_and_scale=True).squeeze().plot()
Out[13]:
<matplotlib.collections.QuadMesh at 0x789d9590da50>
No description has been provided for this image
In [14]:
# Plot the last .tif image
rxr.open_rasterio(ndvi_paths[-1], mask_and_scale=True).squeeze().plot()
Out[14]:
<matplotlib.collections.QuadMesh at 0x789d954794d0>
No description has been provided for this image
In [15]:
# Plot first and last image side-by-side for comparison
recent_ndvi = rxr.open_rasterio(ndvi_paths[-1], mask_and_scale=True).squeeze()
old_ndvi = rxr.open_rasterio(ndvi_paths[0], mask_and_scale=True).squeeze()

# Create side-by-side subplots
fig, axes = plt.subplots(1, 2, figsize=(12,4)) # 1 column and 2 rows

# Plot each in their own axes
old_ndvi.plot(ax=axes[0], cmap=plt.cm.PiYG, vmin=-1, vmax=1) # Set axes between -1 and 1
axes[0].set_title("NDVI - Gila River 2001")
gila_gdf.plot(ax=axes[0], edgecolor='black', facecolor='none',linewidth=1)

recent_ndvi.plot(ax=axes[1], cmap=plt.cm.PiYG, vmin=-1, vmax=1, ) # Set axes between -1 and 1
axes[1].set_title("NDVI - Gila River 2022")
gila_gdf.plot(ax=axes[1], edgecolor='black', facecolor='none',linewidth=1)
Out[15]:
<Axes: title={'center': 'NDVI - Gila River 2022'}, xlabel='x', ylabel='y'>
No description has been provided for this image

Plot NDVI data differences¶

In [16]:
# Compute the difference in NDVI before and after the return of eater rights in 2004
ndvi_diff = (
    ndvi_da
        .sel(date=slice('2012', '2022')) # Choose the dates between 2012-2022
        .mean('date') # Calculate the mean based on the date
        .NDVI 
    - ndvi_da
        .sel(date=slice('2001', '2011')) # Choose the dates between 2001-2011
        .mean('date')
        .NDVI
)

ndvi_diff

#ndvi_diff.plot()
Out[16]:
<xarray.DataArray 'NDVI' (y: 203, x: 382)> Size: 310kB
array([[-0.05567682, -0.0292117 ,  0.00586349, ...,  0.01543377,
         0.01543377,  0.00927271],
       [-0.07940263, -0.03390124, -0.02959213, ...,  0.01815718,
         0.01815718,  0.0177182 ],
       [-0.17723629, -0.08530393,  0.01360923, ...,  0.01517531,
         0.00823637,  0.01195324],
       ...,
       [-0.0115844 , -0.0115844 , -0.00991558, ..., -0.00157142,
        -0.00157142,  0.00205326],
       [-0.01115062, -0.01115062, -0.00994415, ...,  0.00598571,
         0.00598571, -0.00095583],
       [-0.00930774, -0.00849222, -0.01209998, ..., -0.0296714 ,
        -0.02090381, -0.03258313]], shape=(203, 382), dtype=float32)
Coordinates:
    band         int64 8B 1
  * x            (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
  * y            (y) float64 2kB 33.39 33.39 33.38 33.38 ... 32.97 32.97 32.97
    spatial_ref  int64 8B 0
xarray.DataArray
'NDVI'
  • y: 203
  • x: 382
  • -0.05568 -0.02921 0.005863 0.00713 ... -0.02967 -0.0209 -0.03258
    array([[-0.05567682, -0.0292117 ,  0.00586349, ...,  0.01543377,
             0.01543377,  0.00927271],
           [-0.07940263, -0.03390124, -0.02959213, ...,  0.01815718,
             0.01815718,  0.0177182 ],
           [-0.17723629, -0.08530393,  0.01360923, ...,  0.01517531,
             0.00823637,  0.01195324],
           ...,
           [-0.0115844 , -0.0115844 , -0.00991558, ..., -0.00157142,
            -0.00157142,  0.00205326],
           [-0.01115062, -0.01115062, -0.00994415, ...,  0.00598571,
             0.00598571, -0.00095583],
           [-0.00930774, -0.00849222, -0.01209998, ..., -0.0296714 ,
            -0.02090381, -0.03258313]], shape=(203, 382), dtype=float32)
    • band
      ()
      int64
      1
      array(1)
    • x
      (x)
      float64
      -112.3 -112.3 ... -111.5 -111.5
      array([-112.309375, -112.307292, -112.305208, ..., -111.519792, -111.517708,
             -111.515625], shape=(382,))
    • y
      (y)
      float64
      33.39 33.39 33.38 ... 32.97 32.97
      array([33.388542, 33.386458, 33.384375, ..., 32.971875, 32.969792, 32.967708],
            shape=(203,))
    • spatial_ref
      ()
      int64
      0
      crs_wkt :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      inverse_flattening :
      298.257223563
      reference_ellipsoid_name :
      WGS 84
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      geographic_crs_name :
      WGS 84
      horizontal_datum_name :
      World Geodetic System 1984
      grid_mapping_name :
      latitude_longitude
      spatial_ref :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      GeoTransform :
      -112.31041665660531 0.0020833333331466974 0.0 33.38958333034212 0.0 -0.0020833333331466974
      array(0)
    • x
      PandasIndex
      PandasIndex(Index([-112.30937498993873, -112.30729165660559, -112.30520832327244,
             -112.30312498993929, -112.30104165660615,   -112.298958323273,
             -112.29687498993985,  -112.2947916566067, -112.29270832327356,
             -112.29062498994041,
             ...
             -111.53437499000816, -111.53229165667501, -111.53020832334187,
             -111.52812499000872, -111.52604165667557, -111.52395832334243,
             -111.52187499000928, -111.51979165667613, -111.51770832334299,
             -111.51562499000984],
            dtype='float64', name='x', length=382))
    • y
      PandasIndex
      PandasIndex(Index([ 33.38854166367555,   33.3864583303424,  33.38437499700925,
              33.38229166367611,  33.38020833034296,  33.37812499700981,
             33.376041663676666,  33.37395833034352,  33.37187499701037,
             33.369791663677226,
             ...
             32.986458330378234,  32.98437499704509,  32.98229166371194,
             32.980208330378794,  32.97812499704565,   32.9760416637125,
             32.973958330379354,  32.97187499704621,  32.96979166371306,
             32.967708330379914],
            dtype='float64', name='y', length=203))
In [17]:
# Plot the difference
(
    ndvi_diff.hvplot(x='x', y='y', cmap='PiYG', geo=True, 
                     title="NDVI Changes at Gila River Indian Commmunity\n2012-2022 vs. 2001-2011")
    *
    gila_gdf.hvplot(geo=True, fill_color=None, line_color='black')
)
Out[17]:

Conclusion¶

In the figure above, green (a positive score in the differences of NDVI) shows that land got healthier. The pink shows how land got less healthy over time. In the figure we can see that a lot of land in the area outside the Gila River Indian Community got worse over the years, but that has not occurred as much within the GRIC boundary.

These results seem to indicate that when water rights were returned to the indigenous community, there was improvements in the overall health of the land.

Research (https://www4.evergreen.edu/sites/default/files/When_Our_Water_returns_10-25-09.pdf) has shown that when the land was not under indigienous ownership, the river was running dry and there was famine and deprivation within the Indian community. This research aligns with the results we see in the above image.

We can also see in the image, significant dark areas of green where there appear to be agricultural crop growths, indicating imrpovement in local farming in the area.

04: Summarize Results¶

STEP 4: Is the NDVI different within the Gila River Indian Community after the return of water rights?¶

In this section, we compute the mean NDVI inside and outside GRIC boundary.

Try It
  1. Check the variable names - Make sure that the code uses your boundary GeoDataFrame
  2. How could you test if the geometry was modified correctly? Add some code to take a look at the results.
In [18]:
# Compute the area outside the Gila River Indian Community Boundary

outside_gila_gdf = (
                        gpd.GeoDataFrame(geometry=gila_gdf.envelope)
                        .overlay(gila_gdf, how='difference') #Geometry of the Gila River Community's envelope (surrounding area)
)

outside_gila_gdf
Out[18]:
geometry
0 MULTIPOLYGON (((-112.30875 32.96704, -112.3087...
In [19]:
# Check that we collected outside boundary 
outside_gila_gdf.plot()
Out[19]:
<Axes: >
No description has been provided for this image
In [20]:
# Clip data to both inside and outside the boundary
ndvi_inside = ndvi_da.rio.clip(gila_gdf.geometry, from_disk=True)
print(ndvi_inside)

ndvi_outside = ndvi_da.rio.clip(outside_gila_gdf.geometry, from_disk=True)

print(ndvi_outside)
<xarray.Dataset> Size: 47MB
Dimensions:      (x: 379, y: 202, date: 154)
Coordinates:
    band         int64 8B 1
  * x            (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
  * y            (y) float64 2kB 33.39 33.38 33.38 33.38 ... 32.97 32.97 32.97
  * date         (date) datetime64[ns] 1kB 2001-01-14 2001-01-16 ... 2022-01-24
    spatial_ref  int64 8B 0
Data variables:
    NDVI         (date, y, x) float32 47MB nan nan nan nan ... nan nan nan nan
<xarray.Dataset> Size: 48MB
Dimensions:      (x: 380, y: 203, date: 154)
Coordinates:
    band         int64 8B 1
  * x            (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
  * y            (y) float64 2kB 33.39 33.39 33.38 33.38 ... 32.97 32.97 32.97
  * date         (date) datetime64[ns] 1kB 2001-01-14 2001-01-16 ... 2022-01-24
    spatial_ref  int64 8B 0
Data variables:
    NDVI         (date, y, x) float32 48MB 0.6146 0.3796 ... 0.1736 0.2146
In [21]:
# Compute mean annual July NDVI

july_ndvi_inside_df = (ndvi_inside
                       .groupby(ndvi_inside.date.dt.year)
                       .mean(...)
                       .NDVI.to_dataframe())

print('Mean NDVI inside Gila Community \n', july_ndvi_inside_df.head())

july_ndvi_outside_df = (ndvi_outside
                       .groupby(ndvi_outside.date.dt.year)
                       .mean(...)
                       .NDVI.to_dataframe())

print('Mean NDVI outside Gila Community \n', july_ndvi_outside_df.head())
Mean NDVI inside Gila Community 
       band  spatial_ref      NDVI
year                             
2001     1            0  0.199645
2002     1            0  0.177933
2003     1            0  0.187302
2004     1            0  0.176162
2005     1            0  0.238630
Mean NDVI outside Gila Community 
       band  spatial_ref      NDVI
year                             
2001     1            0  0.247629
2002     1            0  0.226726
2003     1            0  0.229889
2004     1            0  0.221753
2005     1            0  0.255275
Mean NDVI outside Gila Community 
       band  spatial_ref      NDVI
year                             
2001     1            0  0.247629
2002     1            0  0.226726
2003     1            0  0.229889
2004     1            0  0.221753
2005     1            0  0.255275
In [22]:
# Join inside and outside data frame 


july_ndvi_df = (july_ndvi_inside_df[['NDVI']]
                .join(july_ndvi_outside_df[['NDVI']],
                      lsuffix=' Inside Gila Community',
                      rsuffix=' Outside Gila Community')) #left is one you start with
print(july_ndvi_df)
      NDVI Inside Gila Community  NDVI Outside Gila Community
year                                                         
2001                    0.199645                     0.247629
2002                    0.177933                     0.226726
2003                    0.187302                     0.229889
2004                    0.176162                     0.221753
2005                    0.238630                     0.255275
2006                    0.211491                     0.235571
2007                    0.181710                     0.211984
2008                    0.201902                     0.237739
2009                    0.179118                     0.218907
2010                    0.201366                     0.235849
2011                    0.173018                     0.217433
2012                    0.180081                     0.222250
2013                    0.181592                     0.223138
2014                    0.192506                     0.224528
2015                    0.200857                     0.227473
2016                    0.179837                     0.213712
2017                    0.195767                     0.225505
2018                    0.178278                     0.210275
2019                    0.212160                     0.227949
2020                    0.205625                     0.224280
2021                    0.223103                     0.237178
2022                    0.197994                     0.212205
In [23]:
# Plot mean NDVI inside and outside the Gila River Indian Community boundary
july_ndvi_df.hvplot(title='Mean July NDVI Inside and Outside Gila River Indian Community')
Out[23]:
In [24]:
 # Add a new column that takes the difference between inside and outside boundaries
july_ndvi_df['Difference'] = (july_ndvi_df['NDVI Outside Gila Community'] 
                              - july_ndvi_df['NDVI Inside Gila Community'])

july_ndvi_df


# Difference is getting less negative, getting closer to zero, getting closer to the same
Out[24]:
NDVI Inside Gila Community NDVI Outside Gila Community Difference
year
2001 0.199645 0.247629 0.047984
2002 0.177933 0.226726 0.048793
2003 0.187302 0.229889 0.042588
2004 0.176162 0.221753 0.045591
2005 0.238630 0.255275 0.016645
2006 0.211491 0.235571 0.024079
2007 0.181710 0.211984 0.030274
2008 0.201902 0.237739 0.035836
2009 0.179118 0.218907 0.039790
2010 0.201366 0.235849 0.034484
2011 0.173018 0.217433 0.044415
2012 0.180081 0.222250 0.042169
2013 0.181592 0.223138 0.041546
2014 0.192506 0.224528 0.032022
2015 0.200857 0.227473 0.026616
2016 0.179837 0.213712 0.033875
2017 0.195767 0.225505 0.029738
2018 0.178278 0.210275 0.031997
2019 0.212160 0.227949 0.015789
2020 0.205625 0.224280 0.018655
2021 0.223103 0.237178 0.014075
2022 0.197994 0.212205 0.014210
In [25]:
# Plot the NDVI difference inside and outside the boundary
july_ndvi_df.Difference.hvplot(title='Difference in NDVI within and outside the Gila River Indian Community')
Out[25]:
In [26]:
%store july_ndvi_df july_ndvi_inside_df july_ndvi_outside_df
Stored 'july_ndvi_df' (DataFrame)
Stored 'july_ndvi_inside_df' (DataFrame)
Stored 'july_ndvi_outside_df' (DataFrame)